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Specialized Monte Carlo simulation techniques and moment free energy method calculations, capable of 
treating fractionation exactly, are deployed to study the crystalline phase behaviour of an assembly of spherical 
particles described by a top-hat "parent" distribution of particle sizes. An increase in either the overall density 
or the degree of polydispersity is shown to generate a succession of phase transitions in which the system 
demixes into an ever greater number of face-centred cubic "daughter" phases. Each of these phases is strongly 
fractionated: it contains a much narrower distribution of particle sizes than is present in the system overall. 
Certain of the demixing transitions are found to be nearly continuous, accompanied by fluctuations in local 
particle size correlated over many lattice spacings. We explore possible factors controlling the stability of the 
phases and the character of the demixing transitions. 



I. INTRODUCTION 

Hard spherical particles can be packed to fill maxi- 
mally just over 74% of space, in the face centred cubic 
(fee) structure. For systems in thermal equilibrium such 
as a colloidal suspension, this structure remains preferred 
^'^ for packing fractions down to about 55%^ where melt- 
ing occurs. But what is the thcrmodynamically optimal 
structure for spherical colloids that are "polydisperse" , 
i.e. have a spread of diameters? Polydispersity should 
act to destabilize a colloidal crystal because of the diffi- 
culty of accommodating a range of particle sizes within 
a single lattice structure; but despite sustained attention 
spanning over three decades (see e.g.^"^'^), there is a lack 
of consensus as to what stable structures arise instead. 

Attempts to address this matter have focused on the 
use of analytical theory and simulation to predict the 
fate of a single crystal in the dense regime (above typi- 
cal fluid densities) when the degree of polydispersity be- 
comes large. Broadly speaking, two incompatible pro- 
posals have emerged: either the system demixes into 
multiple coexisting crystalline phases'"'^ or, alternatively, 
crystalline phases disappear altogether, -"^^ the crystal be- 
ing replaced by an equilibrium glassy phase. ■'^"'^ Ideally, of 
course, one should like to settle the matter as to which 
(if either) of these scenarios is correct by simply perform- 
ing an experiment with a suitable suspension of colloids. 
But the inhibition of diffusion in crystalline phases is ex- 
pected to render solid-solid demixing transitions largely 
unobservablc on experimental timescalcs, even if they are 
thcrmodynamically favoured^**. This would - on the face 
of it - appear to render our central question moot. How- 
ever, one should recognize that even when equilibrium 
is itself imattainable in practical situations, independent 
knowledge of the stable state represents an important 
baseline for interpreting dynamical properties of colloidal 
systems (such as crystallization kinetics^^) which can be 
understood in terms of the topology of the free cincrgy 
surface. There are also suggest icms^'^ that the equilib- 
rium phase diagram sheds light e.g. on the ability of a 
glassy phase to crystallize. The question as to the na- 
ture of the true stable state is thus of more than merely 



academic interest. 

In our view, the disparity in the predictions of pre- 
vious theoretical and computational work is traceable 
to the fact that when considering phase separation, lit- 
tle or no account was taken of "fractionation", i.e. the 
phenomenon whereby the distribution of the particle 
diameters, o", can vary from one coexisting phase to 
another. Indeed it is now well established that frac- 
tionation can radically alter the qualitative features of 
phase behaviour in polydisperse systems compared to 
their monodisperse counterparts (see^" for a review). Ac- 
cordingly, it is essential to fully incorporate its effects if 
one hopes to describe the equilibrium phase behaviour of 
polydisperse systems correctly. 

To quantify fractionation^^ one simply counts, for a 
certain phase (labeled a), the number density of parti- 
cles having diameters in the range a . . .a + da. This 
serves to define a density distribution p^"'^(a). However, 
in real colloidal suspensions, one has the constraint that 
the overall distribution of sizes (across all phases) has a 
form fixed by the synthesis of the suspension. This gives 
rise to a generalized lever rule: 

p«')(a)=$:A(")p(")(a), (1) 

a 

with p''°\a) the "parent" density distribution, p^"^ (a) 
the "daughter" distributions, and A*^"^ the fractional vol- 
ume occupied by phase a (so that A^"^ = 1). Since 
the form of the parent is fixed, only its scale is free 
to vary, e.g. by dilution with solvent, and one writes 
p(o)((T) = n^^^f{a), where n^°^ is the total number den- 
sity and /(cr) is a prescribed normalized shape function. 
The degree of polydispersity, 6, is then defined as the 
standard deviation of the parent distribution /(cr), ex- 
pressed in units of its mean. 

Fractionation greatly complicates the task of determin- 
ing the phase behaviour of polydisperse systems com- 
pared to their monodisperse counterparts. To illustrate 
this, consider (for a given colloidal system) increasing 
j^(o) fj-oni an initially low value, i.e. following a "dilu- 
tion line" through the phase diagram of the system. For 
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sufficiently large n^") the system typically encounters a 
coexistence region of the phase diagram, which is entered 
at a "cloud" point^*^ value of n^'^K At and beyond this 
density the system separates into differently fractionated 
daughter phases. However, as a consequence of fraction- 
ation both the daughter distributions p'"^ (cr) themselves 
and their associated fractional volumes A^") depend non- 
linearly on n^^^. Thus in order to quantify the phase 
behaviour one is faced with the challenge of determin- 
ing the daughter phase properties for all values of n^°^ 
within the coexistence region - a situation which con- 
trasts with the monodisperse case where the coexistence 
densities are independent of the total density, whilst the 
fractional volumes depend linearly on it. 

One theoretical technique that does take fractionation 
into account exactly (within the context of a mean field 
framework) is the moment free energy (MFE) method. 
Previous work using this method by one of us^"^ predicts 
that, for polydisperse spheres, increasing S or n*-"-* within 
the solid region leads to a succession of phase transi- 
tions in which the system demixes into an ever greater 
number of differently-fractionated daughter phases. Each 
daughter phase contains a narrower distribution of par- 
ticle diameters than the parent. This MFE calculation 
thus provides clear evidence for the scenario of multiple 
coexisting solids. But it uses approximate free energy 
expressions, which for solids are derived from those of 
binary mixtures and implicitly already assume that all 
solids are fee. Independent confirmation of its predictions 
is then highly desirable, but has hitherto been lacking. 
The purpose of the present report is therefore to provide 
a definite answer to the question of the nature of the 
equilibrium phase behaviour via state-of-the-art Monte 
Carlo (MC) simulations, and to compare with MFE cal- 
culations; both fully provide for fractionation and employ 
a fixed parent size distribution. 

Our paper is organized as follows. In Sec. II we in- 
troduce our model systems: size disperse hard spheres 
(which we have studied by the MFE method), and 
soft spheres (which we have studied by MC simula- 
tion). Sec. Ill provides a brief description of both the 
MFE method and the bespoke MC techniques required 
for dealing with fixed polydispersity and fractionation. 
Thereafter, in Sec. IV, we report our observations con- 
cerning the phase behaviour of these models, the central 
finding being that the original MFE calculations are in- 
deed correct: as S and/or n*^"-' are increased, a succession 
of transitions occurs in which the system demixes into 
first two, then three, then four fractionated coexisting fee 
crystalline phases. We analyse the observations to arrive 
at a qualitative picture of when a crystalline phase will 
become unstable to phase separation. In Sec. V, we inves- 
tigate in detail the character of these phase transitions, 
finding that some are strongly first order, whilst others 
are quasi-continuous. To quantify the differences, we in- 
troduce and measure a susceptibility that probes parti- 
cle size fluctuations. The associated correlation length at 
the near continuous transitions is found to extend over 



several crystal lattice spacings. Prompted by these ob- 
servations, we use MFE calculations to study in detail 
how the shape of the size distribution affects the ten- 
dency of a given solid phase to exhibit a near continuous 
demixing transition. This leads us to a simple approxi- 
mate criterion that quantifies this tendency. Finally, in 
Sec. VI we summarize and discuss the significance of our 
results, and indicate some issues for future work. 



II. MODELS 

The systems that we shall consider in this work are as- 
semblies of spheres interacting either by a repulsive soft 
sphere potential (as considered by simulation) or a hard 
sphere potential (as studied in our MFE calculations). 
The soft sphere interaction potential between two parti- 
cles i and j with position vectors and rj and diameters 
(7j and aj is given by 

v{rij) = e{GijlrijY^ , (2) 

with particle separation rij = \ri — rj \ and interaction ra- 
dius Gij = {ai + aj)/2. The choice of this potential rather 
than hard spheres is made on pragmatic grounds; in our 
isobaric SGCE simulations (to be reported below), any 
MC contraction of the simulation box that leads to an in- 
finitesimal overlap of two hard spheres would always be 
rejected, so (particularly at high densities) we can expect 
higher MC acceptance rates using a "softer" potential. In 
common with hard spheres, the monodisperse version of 
our model freezes into an fee crystalline structure, ^'^"^^ 
and temperature only plays the role of a scale: the ther- 
modynamic state depends not on and T separately 
but only on the combination n(°)(e/fcBr)^/^. Phase dia- 
grams for different T then scale exactly onto one another, 
and we can fix e/k-QT — 1. 

In all cases we consider parent size distributions of the 
top-hat form: 

..s^ 1(20)-! iil-c<a/a <l + c 

^ \ otherwise ' ^ ^ 

Here the width parameter c controls the degree of poly- 
dispersity 5 = c/\/3. In the following we use the mean 
particle diameter a as our unit of length. 



III. METHODOLOGIES 

A. Analytical calculations: the moment free energy 
method 

Calculating analytically the phase behaviour of poly- 
disperse systems is a challenging problem. This is be- 
cause for each of the infinitely many different particle 
sizes a one has a separate conserved density p{(t). Effec- 
tively one thus has to study the thermodynamics of an 
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infinite mixture, where e.g. from tire Gibbs rule there is 
no upper limit on the number of phases that can occur. 

The moment free energy (MFE) method^^"^° is de- 
signed to get around this issue by effectively projecting 
the infinite mixture problem down to that for a finite 
mixture of "quasi-species" . This is possible when the 
free energy density has a so-called truncatable form, 

f = k^T j dapia) Mpia)) - 1] + ri{pi}) , (4) 

where the excess part depends on a finite number of 
moments of the density distribution, 

Pi = j dup{a)wi{a) . (5) 

This truncatable structure obtains for a large number of 
models of mean field type. Importantly for our purposes, 
it is also found in accurate free energy expressions for 
polydisperse hard spheres, with the simple weight func- 
tions Wi(cr) = (T* (i = 0, 1,2,3). Specifically, we use the 
free energy developed by Bartlett'^^ on the basis of the 
simulation data of Kranendonk et al.^^ for binary mix- 
tures. As mentioned above, this effectively presupposes 
fee structures for all solids, so that validation e.g. by sim- 
ulations, as provided in this paper, is important. 

The MFE method provides a way of expressing the 
ideal contribution to the free energy from Eq. (4) , which 
depends on the complete shape of the density distribu- 
tion, in terms of the moment densities pi. The result is 
the moment free energy. The key feature of the method 
is that if one then treats the quasi-species densities pi as 
if they were densities of ordinary particle species, and 
calculates phase equilibria accordingly, the results for 
cloud points are fully exact. Within coexistence regions, 
the method can be extended by including additional mo- 
ments. Their weight functions can be chosen adaptively, 
and using the resulting approximation as an initializa- 
tion,^^ the exact phase equilibrium conditions can then 
be solved numerically, even if e.g. three or four daugh- 
ter phases are present. Overall, the MFE approach is 
therefore the method of choice for our current investi- 
gation. We do not give further details of the numerical 
implementation here as these are set out in full in Ref..^^ 

B. Simulation: phase behaviour within the isobaric 
semi-grand canonical ensemble 

The appropriate ensemble for determining phase be- 
haviour in dense assemblies of polydisperse particles is 
the isobaric variant of the semi-grand canonical ensemble 
(SGCE).^^ Within this ensemble, the particle number A^, 
pressure p, temperature T, and a distribution of chem- 
ical potential differences fi{a) are all prescribed, while 
the system volume V, the energy, and the form of the in- 
stantaneous density distribution p{(j) all fluctuate. The 
fluctuations in p{u) are linked to the volume fluctuations 



by the relation V ^ p{(T)d(j = N. Importantly, they per- 
mit the sampling of many realizations of the polydisperse 
disorder, thus ameliorating finite-size effects. Moreover, 
in conjunction with volume fluctuations, they facilitate 
separation into differently fractionated phases. Coexis- 
tence of two or more phases is signalled by a multimodal 
form for the distribution of some order parameter such 
as the overall number density n = N/V or the volume 
fraction rj, which for a phase with density distribution 
p{a) can be written as rj = J dap{a){Tr/6)a^ . 

Operationally, the sole difference between the isobaric 
semi-grand canonical ensemble and the constant- A^pT 
ensemble"^^ is that one implements MC updates that se- 
lect a particle at random and attempt to change its diam- 
eter (T to a' by a random amount a' —a drawn from a zero- 
mean uniform distribution. This proposal is accepted or 
rejected with a Mcitropolis probability controlled by the 
change in the internal energy and chemical potential:^^ 

Pace = min [1, exp (-^[A$ + p{a) - /x(<7')])] , 

where A$ is the internal energy change associated with 
the resizing operation and (5 = l/(kBT). 

For SGCE simulations of a polydisperse system at 
some given A'' and T, it is necessary to first determine 
the pressure p and distribution of chemical potential 
differences p{a) such that a suitably defined ensemble- 
averaged density distribution matches the prescribed par- 
ent p'-°^((t) = n^°^f{a). Unfortunately, this task is com- 
plicated by the fact p and jl{a) are unknown junctionals 
of the parent. To solve this problem and hcincc deter- 
mine correct coexistence properties - we shall employ a 
version of a scheme originally proposed in the context of 
grand canonical ensemble studies of polydisperse phase 
coexistence'^^ and later extended to the SGCE,^^''^'' the 
latter implementation of which we now summarize. 

The strategy is as follows. For a given choice of n^°^ 
and temperature T, one tunes p, fl{a) and the A^"^ itera- 
tively within a histogram reweighting (HR) framework,'*'' 
such as to simultaneously satisfy both a generalized lever 
rule and equality of the probabilities of occurrence of the 
phases, i.e. 

n^''^f{<j) = Y^\^'^^p(°'\a), (6a) 

5 = 0, (6b) 

with £ as defined in Eq. (8) below. In the first of 
these constraints, Eq. (6a), the ensemble averaged daugh- 
ter density distributions ^("^(cr) are assigned by aver- 
aging only over conflgurations belonging to the respec- 
tive phase, distinguishable via the multimodal character 
of the order parameter distribution p{n). The devia- 
tion of the weighted sum of the daughter distributions 
p{u) = X)a from the target n^^^f{cr) is con- 
veniently quantifled by a "cost" value: 

A = / I p{a) - n(0)/(c7) | da . (7) 
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In the second constraint, Eq. (6b), 

^-E(.'°'-sy («) 

a ^ ' 

provides a measure of the extent to which the probabihty 
of each phase occuring, , is equal for each of the m 
coexisting phases. Imposing this equality ensures that 
finite-size errors in coexistence parameters are exponen- 
tially small in the system volume. ^^'^"^ 

The iterative determination of p, /i(cr) and A^"^ such as 
to satisfy Eqs. (6a) and (6b) proceeds thus: 

1. Guess initial values of the fractional vohimes A^") 
corresponding to the chosen value of n'^"' . Usually if 
one starts near a cloud point, the fractional volume 
of the incipient phase will be close to zero. 

2. Tune the pressure p (within the HR scheme) such 

as to minimize A. 

3. Similarly tune /i(a) (within the HR scheme) such 
as to minimize A. 

4. Measure the corresponding value of E. 

5. \i£ < tolerance, finish, otherwise vary A(") (within 
the HR scheme) and repeat from step 2. 

In step 3 the minimization of A with respect to varia- 
tions in is most readily achieved'*-^ using the follow- 
ing simple iterative scheme for Ilia): 

/3/ife+i(a)=/3/ife(a) + aln(^-_^j , (9) 

for iteration A; ^ fc -|- 1. This update is applied simulta- 
neously to all entries in the histogram of /i((7), and there- 
after the distribution is shifted so that /i(o'o) = 0, where 
(To is the chosen reference size. The quantity < a < 1 
appearing in Eq. (9) is a damping factor, the value of 
which may be tuned to optimize the rate of convergence. 
Note that (as described in^^) it is important that one 
minimizes A and f to a very high precision in order to 
ensure that the finite-size effects are exponentially small 
in the system size. Typically we iterated until both were 
less than 10~^^. 

The values of A^"^ and p resulting from the application 
of the above procedure are the desired fractional volumes 
and pressure corresponding to the nominated value of 
n^"). As mentioned above, daughter phase densities and 
volume fractions are obtainable by monitoring the mul- 
timodal nature of the order parameter distribution p(n) , 
which allows configurational properties to be assigned to 
a given phase. "^^ 

IV. PHASE DIAGRAM AND SOLID STABILITY 

We consider first the overall phase diagram of the soft 
sphere system as studied in our simulations for a system 



of = 256 particles. Fig. la shows (empty symbols) 
the boundaries of the fluid-solid (FS) coexistence region 
at low densities. These boundaries are the cloud curves 
coming from the low and high density regimes, respec- 
tively, and were previously determined by us using MC 
phase-switch techniques. Oiu- focus in this paper is on 
the solid region at higher densities. Here a comprehensive 
exploration of the (n^^^-S) plane is impractical because of 
the relatively high computational cost of our specialized 
simulation technique. But we can understand important 
qualitative features by following the dashed trajectory 
included in Fig. la'*'^. Along this path, we monitored the 
state of the system via the probability distribution of the 
fluctuating total number density p(n), which serves as 
an order parameter for phase changes. Starting from the 
fee solid cloud point at 5 = 6.3%, we initially increased 
in a stepwise fashion (filled circles) to v!^^^ = 1.45, 
and then switched to increasing 5 at constant n^"^ as a 
potentially faster route to demixing. Indeed, at 5 « 8% 
there was a smooth change in p{n) from single to double 
peaked; an example of the double peaked form is shown 
in Fig. 2a. The two associated phases were identified as 
being fee solids. As is physically reasonable, the higher 
density solid (HDS) daughter phase contains a surplus of 
the smaller particles while the lower density solid (LDS) 
phase has more of the larger particles; see Fig. 3 below. 

Continuing to higher 5 eventually led to spontaneous 
melting of the system at 5 = 13.7%, implying that the 
limit of mctastability with respect to a fluid-solid-solid 
(FSS) coexistence had been overstepped, as is indeed pre- 
dicted by our MFE calculations (see Fig. lb). We there- 
fore backtracked slightly into the solid-solid (SS) region, 
embarking on a new trajectory with increasing rS^^ at 
constant 5 = 13.5%. This produced a third peak in p{n) 
at 71^°^ ~ 1.475. The corresponding intermediate density 
solid (IDS) was again found to be isostructural with the 
other two, with dominant particle sizes between those in 
the HDS and LDS. Finally, increasing the overall density 
to w 1.68 we observed that the central IDS peak 
in p{n) split rather smoothly into two peaks, yielding a 
four peaked structure (Fig. 2b) . All four solids were again 
identified as having an fee structure. 

We next compare to our theoretical MFE calculations. 
These used the same parent size distribution (3) but, 
as explained above, the analysis was performed for hard 
spheres. The reason is that no suitable polydisperse 
model free energies are available for the soft repulsive po- 
tential (2). Nevertheless, the qualitative physics should 
be the same. Indeed, taking a comparable path (see 
Appendix A) through the calculated phase diagram, we 
find the same features as in the simulations, as shown 
in Fig. lb. Quantitatively, the fluid-solid coexistence re- 
gion is narrower, and transitions to multiple solids occur 
at lower n^"^ and S, presumably because with a hard re- 
pulsion, a crystal can accommodate above average-sized 
particles less easily. 

A key feature of the phase diagram is the absence of 
glassy phases. The frustration that could otherwise en- 
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FIG. 1. (Color online), (a) Simulation results for the partial 
phase diagram of the model (2) with parent distribution (3). 
Asterisks: points where new solid phases appear; dashed 
lines: phase boundary slopes found by histogram reweighting. 
F=fluid, S=solid. Colored symbols: state points considered in 
Fig. 3, 4 and 5. (b) MFE calculation of phase diagram of hard 
spheres with the same parent form. The dashed line shows 
a trajectory comparable to that followed by the simulations. 
The cross marks the critical point for the S-SS transition. 



gender such phases is avoided precisely by fractionation. 
To illustrate this, we show in Fig. 3 the density distribu- 
tions for two coexisting solids, at the state points marked 
by the circles in Fig. 1. The figure also shows the par- 
ent density distribution. It is likely that if a single solid 
were forced to have this size distribution at the density 
considered, it would indeed assume a disordered, glassy 
structure. Our results show that at equilibrium, this is 
avoided by effectively splitting the range of particle sizes 
among two phases, allowing each phase to remain crys- 
talline on account of its now narrower range of particle 
size variation. This scenario is then broadly in line with 
that proposed by Bartlett,^ but the split in sizes is not 
"sharp" in the sense that particles of a given size would be 
found exclusively in one phase or the other. Such a sharp 
split would require infinite differences between phases of 
the relevant size-dependent chemical potentials. Apart 
from the general phenomenon of fractionation. Fig. 3 
also demonstrates good agreement between the simula- 
tion results and the MFE predictions, with e.g. the cross- 



FIG. 2. Order parameter distributions through (a, top) the 
S-SS transition and (b, bottom) the SSS-SSSS transition. 



ing point between the three density distributions in both 
cases located somewhat to the right of the parental mean. 

As the density or parent polydispersity of a system in 
the SS regime are increased further, the polydispersity in 
the two daughter phases becomes unfavorably large. At 
this point a third solid appears that takes up the middle 
of the size distribution, producing three daughters whose 
size distributions are again sufficiently narrow. This is 
illustrated in Fig. 4. At the transition from this SSS 
regime to four solids (SSSS), we then see a process that 
is qualitatively similar to the S-SS transition: the middle 
(IDS) phase splits into two phases, each again with a nar- 
rower size distribution (Fig. 5). It is worth emphasizing 
that also in these more complicated fractionation scenar- 
ios, the agreement between simulations for soft spheres 
and theory for hard spheres remains good. 

A natural question to ask about the results so far is: 
what determines the stability of solid phases, i.e. when 
do new solids appear? Intuitively one would expect that 
there should be a certain threshold in polydispersity be- 
yond which a given single solid phase would become 
thermodynamically unfavourable. This threshold should 
then depend on how dense the phase is: a denser solid 
can accommodate less variation in particle sizes. To test 
this idea quantitatively, we plot along the path through 
our phase diagram the polydispersity 6 versus the volume 
fraction rj of all coexisting phases. The results are shown 
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FIG. 3. Density distributions in the SS regime, (a) Solid line: 
Parent density distribution at state point (n'°' = 1.45,5 = 
9.5%), marked by the red circle in Fig. la. Symbols: Simu- 
lation results for the two daughter distributions. The associ- 
ated fractional volumes A'"' are 0.527 (HDS, squares), 0.473 
(LDS, circles), (b) MFE results at the comparable state point 
(n(°^ = 1.133, 5 = 6.12%), marked by the red circle in Fig. lb. 
Fractional volumes are 0.561 (HDS), 0.439 (LDS). 



FIG. 4. Density distributions in the SSS regime, (a) Solid 
line: Parent density distribution at state point (n^"' = 
1.60,(5 = 13.5%), marked by the green triangle in Fig. la. 
Symbols: Simulation results for the three daughter distribu- 
tions. The associated fractional volumes A'"-* are 0.267 (HDS, 
squares), 0.309 (IDS, diamonds), 0.424 (LDS, circles), (b) 
MFE results at the comparable state point (n^"' = 1.186, 
5 = 8.7%), marked by the green triangle in Fig. lb. Frac- 
tional volumes are 0.341 (HDS), 0.253 (IDS), 0.405 (LDS). 



in Fig. 6. One sees that the coexisting phases do indeed 
cluster around a line in the (77, S) plane, although the 
clustering is clearly tighter for the MFE (hard sphere) 
theory. In the plot for the latter case we also show the 
S-SS phase boundary from Fig. lb as a dashed line. Re- 
call that this is the boundary as it applies to solids with 
a top hat size distribution. Most of the coexisting phases 
that we find lie inside this phase boundary, implying that 
with their smoother size distributions they can tolerate a 
somewhat larger amount of polydispersity. In summary, 
while the general picture of a line in the volume fraction- 
polydispersity plane where solids become unstable holds 
true, this line is broadened into a transition region by its 
dependence on the shape of the size distribution. Mo- 
tivated by this finding we also experimented with other 
measures of polydispersity to see whether they would re- 
duce this dependence on distribution shape. In particu- 
lar, we considered d2n = [^((o- - a')^")]^/(^"V('7> where 
the averages are over particle sizes a and a' randomly 
drawn from the relevant size distribution. For n = 1 
this gives the conventional S: for n ^ 00 it becomes the 
difference between the largest and the smallest particle 
size present, normalized by the mean size. While one 
may imagine the latter quantity to be the most relevant 
one for determining crystal stability, we found in practice 
that the clustering in the (rj, d2n) plane becomes worse for 
larger n, with the most easily interpretable results being 
the ones shown above for n = 1. 



As a final comment on Fig. 6 it is worth highlighting 

that the phase with the highest density (HDS, shown by 
squares) in fact always has the smallest volume fraction 
among the daughter phases for a given parent. This again 
reflects the strong fractionation effects: as illustrated in 
Figs. 3-5, the HDS phase contains the smallest particles, 
and this reduces its volume fraction 77 to the point where 
it is smaller than for all other phases. This trend is true 
throughout, i.e. the ordering of the daughter phases by 
density n is always the reverse of the ordering by volume 
fraction 77. 

V. CRITICALITY IN TRANSITIONS TO MULTIPLE 
SOLIDS 

A. Order parameter distributions and fractional volumes 

In this section we discuss the nature of the transitions 
as our system of polydisperse spheres fractionates into 
an increasing number of solids. Our focus will be on the 
rather surprising finding that these transitions can be 
nearly continuous in character. 

Initial evidence for this claim is provided by Fig. 2 
above. This shows the distributions p{n) of the fluctu- 
ating number density in the MC simulations, with each 
peak corresponding to one of the solid phases. One sees 
in Fig. 2a, for the S-SS transition, that the initial sin- 
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FIG. 5. Density distributions in the SSSS regime. (a) 
Solid line: Parent density distribution at state point (n'"^ = 
1.73, 5 = 13.5%), marked by the blue square in Fig. la. Sym- 
bols: Sirrmlation results for the four daughter distributions. 
The associated fractional volumes A'"' are, from left to right: 
0.209 (HDS, squares), 0.188 (IDS2, +), 0.232 (IDSl, dia- 
monds), 0.373 (LDS, circles), (b) MFE results at the compa- 
rable state point (n'"-' = 1.186, 5 = 8.7%), marked by the blue 
square in Fig. lb. Fractional volumes are, from left to right, 
0.273, 0.162, 0.200, 0.365. From Ref..*'* Copyright American 
Physical Society. 

gle peak splits smoothly into two nearby peaks which 
then rapidly move outwards towards more clearly sep- 
arated densities. This contrasts with what one would 
have expected for a first order transition, where a new 
peak appears at some finite distance from the initial peak 
and gradually acquires more and more weight. Such 
a scenario is found, along our particular path through 
the phase diagram, for the SS SSS transition (data not 
shown). The SSS-SSSS transition, on the other hand, is 
again nearly continuous, like the S-SS transition. This 
can be seen in Fig. 2b, where the middle peak splits 
smoothly into two new peaks which move apart and form 
the IDSl and IDS2 phases. 

Further evidence for nearly continuous transitions to 
multiple solids is provided by the variation of the frac- 
tional phase volumes X^°\ shown in Fig. 7. One observes 
that at the S-SS transition, the fractional volume occu- 
pied by the new phase has a strongly nonlinear variation 
with the parent polydispcrsity. In fact, looking at the 
simulation results (Fig. 7a), where we cannot get reli- 
able data close to the transition, one would guess that 
the fractional volume of the new phase has a discontin- 
uous onset, as is typical of phase transitions which are 
continuous in the thermodynamic sense. The difficulty 
in obtaining data close to the transition in simulations 
stems from the fact that in a finite-sized system the crit- 
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FIG. 6. (a, top) Location of daughter phases along the verti- 
cal and final horizontal paths in the phase diagram of Fig. la, 
plotted in terms of volume fraction rj and polydispersity S. 
Solid, gray and empty symbols refer to the SS, SSS and SSSS 
regions, respectively, (b, bottom) Analogous plot for MFE 
calculations for hard spheres. The additional dashed line indi- 
cates the S-SS cloud curve for top hat size distributions from 
Fig. lb. 



ical density distribution p(n) has two peaks, so that one 
has to proceed some way into the two phase region be- 
fore one can be sure that peaks observed in p{n) indicate 
genuine phase coexistence. Looking at the right half of 
Fig. 7, the behaviour at the SS-SSS transition is rather 
different, with the fractional volume taken up by the new 
phase increasing smoothly from zero in an almost linear 
fashion. This is in line with expectations for a first or- 
der transition. The SSS SSSS transition, on the other 
hand, again shows nearly continuous behaviour. As for 
the S-SS transition, phase coexistence cannot be deter- 
mined unambiguously from the simulation data for our 
finite systems, and the data outside of the resulting gap 
are again suggestive of a jump in the new fractional vol- 
ume at the transition. The MFE calculations show that 
there is no real jump, rather a strongly nonlinear increase 
from zero, so that the transition is close to but not fully 
critical. 

Taken together, the above observations of the be- 
haviour of the density distribution p{n) in the simula- 
tions, and of the variation of the fractional phase vol- 
umes, provide strong evidence that demixing transitions 
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FIG. 7. (a, top) Variations of A through transitions from 
single soHd to multiple solids (SS, SSS, SSSS). (b, bottom) 
Corresponding results from MFE calculations. 



to multiple solids can be near critical. Along our specific 
path through the phase diagram, it is the S-SS and SSS- 
SSSS transitions that are of this type. To investigate this 
issue in more detail, we now turn to characterizing the 
near critical properties at the level of single phases, via 
appropriate correlation functions of particle size fluctua- 
tions. 



B. Correlations in size fluctuations 

To define a measure of how strongly spatially corre- 
lated size fluctuations are in our solids, we consider first 
a grand canonical setting for a single phase in a fixed 
volume V, and with imposed chemical potentials i^{a). 

The fluctuating density distribution p{a) has ensem- 
ble average (pier))- If we define moment densities p^i — 
/ dap{a)a", then the normalized ensemble average size 
distribution is {p{a))/{po). Its variance E = {p2)/{po) — 
((pi)/(po))^ sets the scale for any particle size fluctu- 
ations. To define our correlation measure x, we mea- 
sure the mean particle size in any configuration, which 
is pi/po, and construct its variance across the ensemble. 
This is then normalized by E and multiplied by system 



volume V to get a quantity with the dimension of a vol- 
ume: 



V{[A{p,/po)]') 



(10) 



In the thermodynamic limit of large V, po and pi 
have small fluctuations so one can expand A(pi//9o) = 
(Api)/(/9o) - (Apo)(/3i)/(po)^- Abbreviating the 
ensemble- averaged mean size as a = {pi)/{po), this gives 



X = 



V{[Ap, - aApoY) 



(po)^S 

or in terms of the fluctuating density distribution 

V{{Jdaia-a)Ap{a)Y) 



X = 



(11) 



(12) 



The denominator here could also be written as {p2) {po) — 
{Pi?- 

To motivate further the above definition of our measure 
of correlations x, one can express it via correlation func- 
tions of the full spatially-resolved density p{r,a). The 
fluctuations of the latter can be expressed in terms of 
the pair correlation function g^^r' {r) between particles of 
sizes a and a' as^^ 

{Ap{r,a)Ap{r',(j')) = {p{(T))S{r' - r)5{a' - a) 

+ (p(c7))(p(a'))[w(r-'-r)-l] 

So the numerator of (12) is, using Ap(cr) = 
V-^J drAp{r,a), 

V-^ j drdr'dadc7'{(7-a){a'-a){Ap{r,a)Ap{r',a')) = 

(po)S+ / drd<jda\a-<j){a'-a){p(a)){p{a'))[g„^,{r)-l] 



This shows that our deflnition of x is physically reason- 
able: it is the volume integral of a correlation function 
that measures the spatial correlations of fluctuations in 
particle size away from the ensemble mean. We will 
therefore also refer to x as the size fluctuation suscep- 
tibility. Note that the trivial first term above makes 
a contribution of l/(po) to x- This is the unit volume 
per particle and of order unity in the density range we 
are considering. Wc will see below that it is negligible 
compared to the main contribution from the correlation 
function integral. 

Some care is needed when relating the susceptibility x 
as defined above to a length scale ^ for the spatial cor- 
relations of size fluctuations. Away from criticality, and 
in d spatial dimensions, then since the correlation func- 
tion being intcigrated decays on a spatial scale of ^, one 
estimates X ^ C^- This is the identification we made 
previously.^'' At criticality, on the other hand, the corre- 
lation function appearing above will have a spatial power 
law decay with |r|~''+^~'' up to the cutoff, and hence the 
susceptibility scales as x ~ where r] is the standard 
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critical exponent (and not, as elsewhere in the paper, the 
volume fraction). 

One can show that, for large systems, the size fluctu- 
ations we are considering are the same in all reasonable 
ensembles, for example a semi-grand canonical ensemble 
where particle number N is fixed and the volume V can 
fluctuate. In this case the factor V in (10) is replaced by 
^/ (po) to give 



X ■ 



N{[A(p,/po)f) 
{P2) - (pi)V(Po) 



(13) 



and this is the method we use to extract x from simu- 
lation data. For the theoretical calculations, we employ 
(11) and extract the (co-) variances of the fluctuations of 
the moment densities po, pi from the appropriate curva- 
ture matrix of the moment free energy. ^'^ 

Note flnally that in the context of experiments on col- 
loids a canonical ensemble, with fixed particle number 
N, volume V and parent size distribution, would be the 
most natural description. For a single phase, the mean 
size is then fixed and no size fluctuations occur. But 
X can still be defined in terms of the pair correlation 
function gacr'{r) as described above, provided the spa- 
tial integration over r is cut off at some distance much 
larger than the correlation length but much smaller than 
the system size. This eliminates the contribution from 
the nonzero values g„„i{r) — 1 = 0{1/N) that remain 
at larger r when the total particle number N is fixed. 
Once several phases appear, each phase has fluctuating 
particle numbers and volume, but one can check that the 
size fluctuations in each phase, and hence the size fluctu- 
ation susceptibility, are as would be calculated for single 
phases in the grand canonical ensemble. 

Having defined how we will quantify the strength of 
correlations in spatial particle size fluctuations, we show 
in Fig. 8 results for x along the vertical and flnal hori- 
zontal paths through the phase diagrams of Fig. 1. One 
observes that x grows large near the transitions to two 
and four solids, confirming their near continuous charac- 
ter. In the latter case, the splitting of the middle peak 
seen earlier in p{n) suggests that the new solids arise out 
of the IDS phase, and this is consistent with large fluctu- 
ations occurring (see Fig. 8) only in this phase and not 
the HDS or LDS. The MFE predictions are, again, in 
good qualitative accord with the simulation data. 

To summarize our observations in this section, the be- 
haviour of fractional phase volumes and of the order pa- 
rameter distributions p{n) suggested that phase transi- 
tions to multiple solid phases can be near critical in na- 
ture; for our path throiigh the phase diagram this applies 
to the S-SS and SSS-SSSS transitions. We proposed the 
size fluctuation susceptibility x as a quantitative measure 
of the range of correlations in the spatial fluctuations of 
particle sizes. Results for this from both simulations and 
MFE calculations then demonstrated that these transi- 
tions are indeed close to critical, being characterized by 
values of x far above the unit volume per particle. 
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FIG. 8. Size fluctuation susceptibility x in the solid phases 
encountered along the phase diagram trajectories of Fig. 1. 
(a) Simulations, (b) MFE calculations. EYom Ref..** Copy- 
right American Physical Society. 



That such critical or near critical transitions from one 
to several solids might occur is plausible given that S-SS 
critical points arc observed also in simulations of binary 
hard sphere mixtures. The free energy expression for 
polydisperse hard spheres that we use in our MFE cal- 
culations was devised by Bartlctt'^^ on the basis of free 
energies fitted to these binary mixture simulations. The 
polydisperse system must then "inherit" the existence of 
critical points, though not in any trivial way. For ex- 
ample, the polydisperse system has many more degrees 
of freedom for fiuctuations in its size distribution, and 
one can show from this that spinodal densities are al- 
ways lower in the polydisperse than in the corresponding 
binary case. 

It is worth stressing that even in transitions involving 
multiple solid phases (SSS-SSSS), criticality is essentially 
a single-phase property. Indeed, we have found above 
quite distinct values of x i^i the three coexisting solid 
phases before the transition to four solids. In the simula- 
tions, we have only a single phase in the simulation box 
for most of the time, which emphasizes further that x is 
determined from the properties of this single phase. To 
be more precise, there is an effect of the presence of other 
phases: the lever rule forces the density distributions of 
all phases to add up to the parent, and this provides con- 
straints on the chemical potentials ij{cr). Once we know 
these chemical potentials, however, we can determine x 
individually for every phase, independently of the others. 

We next ask what features of a given size distribution 
make it undergo a critical or near critical transition to 
multiple solids. Having recognized that criticality is a 
single-phase property, we focus in this enterprise on the 
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cloud point for the transition S-SS from a single solid to 
two fractionated solids. 



C. Predicting criticality 

The question of determining whether the S SS transi- 
tion from a parent with a given size distribution is close to 
critical can be cast in quantitative terms as follows: how 
large is the size fluctuation susceptibility x s-t the S-SS 
cloud point? We investigate this using MFE calculations 
for the hard sphere case; precise simulation studies would 
inevitably require finite-size scaling to larger system sizes 
than we can access using our computational resources. 

As explained in Sec. Ill A, the free energy expression 
that we use for polydisperse hard spheres have excess 
contributions that depend only on the moment densi- 
ties Pi with I = 0,1,2,3, defined by the weight func- 
tions Wi{a) = cr*. One can then show in generality (see 
e.g.^^) that the criterion for a spinodal, where a phase 
becomes unstable to local density fiuctuations, involves 
these moments as well as those defined by the second- 
order weight functions Wi{a)wj{a) = f7'+-', giving in our 
case moments up to pg. For a given particle size distribu- 
tion, all the ratios pi/po, . . . jPq/pq are fixed and the den- 
sity po = n at the spinodal can be found from the spin- 
odal criterion. The additional condition for a spinodal 
point also to be a critical point involves in addition the 
third-order weight functions iVi(a)wj{(j)wk{(y) = (T*^-'^'^, 
which produce moments up to pg. Inserting the spin- 
odal density, the exact critical point condition resulting 
from our model free energies is then some function of 
pi/po, ■ ■ ■ , pg/po- These are the 1st to 9th moments of 
the normalized size distribution, and so whether a parent 
phase with a given size distribution will exhibit a critical 
S SS transition or not depends only on these moments. 

Unfortunately, because the solid free energies we use 
are derived from fits to simulation data,'^^ the critical 
point condition that results is far too complicated to al- 
low for any analytical progress. We therefore proceed 
initially by solving the condition numerically for a range 
of parent size distributions of interest. The first case 
to consider is evidently the top hat distribution studied 
throughout the paper so far. Here there is only a single 
parameter to vary, namely the polydispersity S. We solve 
for each S the spinodal condition to find the spinodal den- 
sity, and then evaluate the critical point condition at this 
density. It turns out that there is indeed a critical point 
in the phase diagram, at (n(°) = 1.1669,(5 = 0.0472). It 
is marked in Fig. lb, and lies close to the path through 
the phase diagram that we have considered above. This 
rationalizes why the S-SS transition along this path is 
near critical, with a large value of x- at the critical point 
itself, we would have found x diverging to infinity at the 
transition. 

The situation with regard to the shape of the parent 

distribution is not trivial, however. For example, in pre- 
vious work we considered both triangular and Schulz dis- 



tributions, and found no critical points on the S-SS 
cloud curve in the physically relevant ranges of density 
and polydispersity. To get more insight, we consider next 
families of parent distributions where we can tune both 
the width, as measured by 5, and the shape. Generaliz- 
ing from the top hat case studied above, we look first at 
"slanted top hat" parents where the size distribution is 
/(ct) = A + Ba in some interval cr_ < <J < <j+, and zero 
otherwise. We adjust A, S, (t_ and so that /(cr) is 
normalized, has mean 1 as before, and the desired value 
of 5. This leaves one degree of freedom, which we express 
via the slant ratio R = /((t+)//((t_), with R = 1 giving 
back the simple top hat distribution. 




FIG. 9. Critical polydispersity 5 versus slant ratio R for 
slanted parents. Dashed: approximation from Sv = 0. 

Proceeding as for the top hat parent, we can now deter- 
mine numerically for fixed slant ratio R the critical value 
(if any) of S, or vice versa. In the resulting Fig. 9 we 
observe that whether or not there are critical points for 
a given parent shape depends on R: for R below around 
0.81, no critical points appear; for slightly larger values, 
two critical points can exist in the phase diagram, and 
for values of R around unity and above we generically 
find one critical point. 

So far we have asked what marks out parent size dis- 
tributions that have critical S-SS transitions, which cor- 
responds to X = oo at the cloud point. Here we digress 
slightly to ask how x at the cloud point then varies as 
we move away from the critical parent shape. Data from 
MFE calculations are shown for this in Fig. 10, where we 
consider parents with fixed density n^^^ = 1.133 as on the 
vertical path in Fig. lb. For given slant ratio R we find 
the polydispersity 6 at the cloud point; see the bottom 
inset of Fig. 10. The main plot displays the resulting 
cloud point value of the susceptibility x against R. It is 
seen to diverge as a critical value R = Rc is approached, 
and indeed by solving the critical point criterion for the 
given parent density we find a single such critical value, 
Rc = 1.734. This means that if we had considered a par- 
ent with this slanted shape, we would have seen - within 
our MFE calculations for hard spheres - a fully critical 
S-SS transition on the vertical path in Fig. lb. 
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FIG. 10. Size fluctuation susceptibility x at S-SS cloud point 
for slanted parents with density n'-'^-' = 1.133, versus slant 
ratio R. Bottom inset: 5 at the cloud point. Top inset: Semi- 
log plot of X vs the deviation from the critical slant ratio, 
\R — Rc\- The dashed line is a power law with exponent —2. 



The top inset of Fig. 10 plots the susceptibility x ver- 
sus the distance from the critical parent shape. The data 
are consistent with a divergence as x ~ l-R ~ -Rc|~^) 
except for the points nearest Rc either side, where our 
numerics become unreUable. That exponent value may 
seem surprising at first: for our mean field free en- 
ergy, the susceptibility for models in the Ising univer- 
sality diverges as x ^ |T — Tc|~^ with 7=1. But 
R smoothly changes the parent shape, and the latter is 
analogous to the Ising magnetization m. We should then 
identify \R — Rc\ with m and this leads to the scaling 
X ~ m~'^/^ ~ \R — Rc\~^^^ ■ For our mean field free 
energy this gives an exponent value 7//? = 1/(1/2) = 2, 
exactly as observed. A more accurate theory which cap- 
tures the non-mean field Ising singularities would then be 
expected to give in d = 3 the exponent 7//? = 5—1 « 3.9. 

The bottom inset of Fig. 10 shows the value of the 
polydispersity 5 at the cloud point against the slant ratio 
R, for the same fixed parent density as in the main plot. 
The variation in 5 is very small, between around 0.052 
and 0.056, even though the parent shape changes quite 
dramatically from i? = 1/4 to i? = 4. This is in line 
with the expectation that 6 is the main aspect of the size 
distribution that determines solid stability. 

Returning now to the question of what determines 
whether a given particle size distribution will produce 
a critical S SS transition, we broaden our investigation 
to a wider class of distributions, namely the Beta distri- 
butions. These are of the form /(cr) oc {a — a-Y{(j+—(j)^ 
in some interval a- < cr < (t+, and zero otherwise. The 
values of the smallest and largest sizes <j- and cr+ and 
the proportionality coefficient are again adjusted to make 
/((j) normalized with unit mean and standard deviation 
5. The advantage of Beta distributions is that with their 
two shape parameters, a and h, they are more flexible 
than e.g. the slanted top hat parents from above. In 
particular, they can interpolate from distributions with 
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FIG. 11. Critical line in the (o, &)-plane for Beta size distribu- 
tions. The density is fixed to n = 1.133 and the polydispersity 
5 is determined from the spinodal condition; it varies very lit- 
tle over the range shown, from 0.0559 to 0.0576. The line 

resulting from the criterion s^, = is also shown. Both lie sig- 
nificantly below the line for symmetric distributions (a = 6). 

fairly sharp cutoffs at the extreme sizes - for low a and b, 
where in particular a = & = gives back a top hat distri- 
bution - to ones with almost Gaussian shape (large a and 
h) where the cutoffs are in the far tails of the distribution. 

We can now proceed as above and solve the MFE crit- 
ical point criterion to find out what shape parameters a 
and h produce critical S-SS transitions. We do this at 
fixed density, taking again n^") = 1.133. The spinodal 
condition then fixes 5 for given a and 6, and the criti- 
cal point gives one additional condition, so that we get a 
line of critical points in the a, h plane as shown in Fig. 11. 
What is noticeable is that the critical points lie signifi- 
cantly away from the line a = h where the; parent density 
distribution is symmetric. This is also clearly visible in 
Fig. 12a, with the critical size distributions having dis- 
tinct peaks to the right of the mean. One is led to ask 
whether there are other quantities, related but not identi- 
cal to particle diameter, that would have more symmetric 
distributions. An obvious choice is the particle volume, 
which is proportional to 1; = cr'^. As Fig. 12b shows, the 
distributions /„(w) = /((7)/(3(t^) are indeed much more 
nearly symmetric at criticality. This suggests that devi- 
ations from such symmetry, as measured by the skew 



indicate deviations from criticality, and conversely Sy = Q 
might be a reasonable approximate way of identifying 
critical size distributions. We have included the line in 
the (a, h) plane that results when we solve this condition 
(at the same values of 5 as previously) in Fig. 11. The 
agreement with the critical line calculated directly from 
the MFE criticality condition is qualitatively quite good. 
In particular, the criterion s„ = captures the fact that 
the critical size distributions are asymmetric when ex- 



12 



a = 0.5 




s ip 

V 

0.8- 
0.6- 
0.4 - 
0.2 - 
0- 
-0.2- 
-0.4- 
-0.6- 
-0.8- 



■ 


s 


□ 


HDS 


o 


LDS 





0.8- 
0.6- 
0.4- 
0.2- 
0- 
-0.2- 
-0.4- 
-0.6- 
-0.8, 



(») 



□□□□□□□□ 



■ 


s 


□ 


HDS 





LDS 



□ 


HDS 





IDS 


o 


LDS 



04 0.05 0.06 0.07 g 0.0 



1.15 (0) 1.2 



FIG. 12. (a, top) Examples of critical Beta size distribu- 
tions, corresponding to the values of (a, 6) from Fig. 11, with 
a = 0.5, 1, 1.5, . . . , 5 increasing in the direction shown, (b, 
bottom) Corresponding distributions of w = cr^, which is 
proportional to particle volume; these distributions are much 
more nearly symmetric. 



FIG. 13. Skewness s„ of particle volume distribution in 
the phases occurring along the phase diagram trajectories of 
Fig. 1. (a, top) Simulation results, (b, bottom) MFE cal- 
culations. Comparison with Fig. 8 shows that the near critical 
phases also have small Sy. 



pressed in terms of particle size, with a > b throughout. 

We have also calculated s„ along the path through the 
phase diagram in Fig. 1 for top hat parents, and show 
the results in Fig. 13. One observes that the parent phase 
has relatively low s„ at the S-SS transition, in agreement 
with the large vahies of the size fluctuation susceptibility 
X- Likewise, in the SSS-SSSS transition, the phase that 
exhibits large x and splits in a near critical fashion into 
two solids also has small Sy. 

Further support for the use of s„ = as an approxi- 
mate criterion for criticality comes from the fact that Sy 
can be written in terms of moment densities of tr as 



{v'^)-3{v^){v) + 2{vy 

- (t;)2)3/2 

3p6P3Po + 2pi 



P9Po 



{pePo - Plf^ 



(15) 
(16) 



which entails exactly the moment densities po,...,pg 
(though not all of them) that we would expect from 
the general discussion above. Nevertheless the criterion 
Sy = clearly remains approximate: for the slanted top 
hat parents, the results in Fig. 9 show that here the agree- 



ment with the full criticality criterion is less good. In 
particular, from s„ = we would predict that there are 
no critical points for slant ratio i? < 1, whereas in fact 
critical size distributions exist down to R 0.81. The 
question of whether there is a more accurate yet still sim- 
ple criterion for S-SS criticality remains open. 



VI. DISCUSSION AND FUTURE WORK 

In summary we have deployed tailored Monte Carlo 
simulation methods and moment free energy calculations 
to provide conclusive evidence that dense polydisperse 
spheres at equilibrium demix into coexisting fee phases, 
with more phases appearing as the spread of diameters 
and the number density increase. Up to four coexisting 
phase were tracked, each of which contained a narrower 
distribution of particle sizes than is present in the system 
overall. Interestingly it was observed that for our systems 
the S-SS and the SSS-SSSS transitions are quasi-critical, 
characterised by a large correlation length for fluctua- 
tions in local particle size. By contrast the SS-SSS tran- 
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sition was found to be strongly first order. To rationalize 
these observations, we investigated the features of the 
parental size distribution that control the eharaeter of 
solid demixing transitions. It was found that small skew 
in the parent distribution of particle volumes (s„ « 0) 
correlates well with the existence of a quasi-continuous 
transition, at least for one class of parental distribution 
shapes. 

Whilst our results settle the matter of the true equi- 
librium behaviour, they leave open the question as to the 
extent to which this behaviour will be observable in ex- 
perimental studies of polydisperse systems. Initial indi- 
cations from recent experiments on colloid-polymer mix- 
tures are that solid-solid demixing does not occur on the 
timescale of weeks. Thus the best opportunity to see ev- 
idence may be to focus on regions of the phase diagrams 
where polydisperse solid(s) coexist with a fluid that can 
transport particles to their preferred solid phase. Addi- 
tionally it would be interesting to try to manufacture a 
distribution of particle sizes that has s„ « and then 
look for an increase in particle size fluctuations in the 
single solid region, even if the full transition itself is not 
seen. 

As regards the questions that our results pose for fur- 
ther simulation and theoretical work, an interesting mat- 
ter is that of the fate of the regions of multiple solid co- 
existence at high volume fraction. As the polydispersity 
S is reduced, it seems clear that all the transition lines 
to multiple solids (S-SS, SS-SSS, SSS-SSSS etc) must 
converge on (but never quite reach) the monodisperse 
close packed limit at ^ = 0, 77 = tt/vTs rs 74% since the 
close packed crystal will be unstable to any finite dc^grec 
of polydispersity. We can also consider what happens if 
we fix the polydispersity 5 > and increase the parent 
volume fraction. The number of fractionated solids will 
increase without bound as the pressure increases, until 
at some volume fraction the pressure diverges and the 
system cannot be compressed further. The locus of these 
points in the phase diagram forms the infinite pressure 
line. Also this line must, as S is decreased to zero, ap- 
proach the monodisperse close packed limit rj 74%. 

An intriguing question is whether criticality can play 
a role in the approach to the close packed limit along 
the S-SS boundary. As we have seen, it is primarily 
the parent shape that controls the nature of the S-SS 
transition. Thus there may exist parent forms for which 
S-SS demixing is critical at or very near to the close 
packed limit, and it would be interesting to see whether 
a simple characterization of such narrow critical parent 
size distribution forms can be found. 

Further solid phases are likely to arise in the phase 
diagram at values of S beyond those that we have ex- 
plored. For example one could imagine that at very large 
S (for which the system separates into multiple coexisting 
phases) the smallest particles, rather than forming their 
own fee phase, might instead secrete themselves in the 
interstitials of the fee solid formed by the largest parti- 
cles, thus potentially permitting the volume fraction to 



exceed 74%. Indeed this could be a mechanism whereby 
a unimodal parental distribution might produce familiar 
substitutionally-ordered phases, such as CsCl, or exotic 
phases such as AB2 and AB13 that can appear in binary 
colloid mixtures.'** Investigating this question could - in 
principle - be tackled by simulation, but is probably be- 
yond the present capabilities of the MFE calculations 
which are based on free energies that are reliable only for 
small to moderate S. 
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Appendix A: Comparable locations in soft and hard sphere 
phase diagrams 

It is difficult to map from first principles the simulation 
phase diagram for soft spheres to the MFE calculations 
for hard spheres. Existing approaches as summarized in 
e.g. Ref.^^ do allow one to calculate effective hard sphere 
diameters for soft particles, but are based on liquid-state 
correlations and work only up to moderate densities. 

We therefore identified comparable points based on 
the phase diagram topology. In particular, the simula- 
tions at polydispersity d = 13.7% show an instability 
towards FSS coexistence very close to the SS-SSS transi- 
tion (see Fig. la). From the density range between these 
two points, relative to the separation between the SS-SSS 
and SSS-SSSS transitions, we estimate the correspond- 
ing polydispersity for hard spheres to be 5 = 8.7%, just 
below the meeting point of the SS-FSS and SS -SSS lines 
in Fig. lb. We find the density corresponding to the 
vertical trajectory in Fig. la (n^°-* = 1.45) similarly: at 
this density and at S = 13.7%), the simulations show an 
SS phase split that is still stable but becomes unstable 
at slightly higher i5. The corresponding density in the 
MFE phase can be estimated as n^^^ = 1.133, just be- 
low the SS-FSS transition line a.t 5 = 8.7%. This fixes 
the vertical and final horizontal trajectories through the 
MFE phase diagram which we use in evaluating e.g. the 
correlation volume data in Fig. 8. 

For the SSSS state point in Fig. 5 we proceed sim- 
ilarly. This point lies on the final horizontal trajec- 
tory through the phase diagram, for which we already 
have the hard sphere polydispersity 6 = 8.7% that cor- 
responds to the simulation value S = 13.7%. We then 
estimate the density of the state point so that its den- 
sity difii'erence to the SSS-SSSS transition, in units of 
the separation between the SS-SSS and SSS-SSSS tran- 
sitions, is the same as in the simulations. This gave 
= 1.232. The same method was applied for the SSS 
state point in Fig. 4. For the SS point in Fig. 3, we sim- 
ply scaled the polydispersities in proportion to the value 
of 6 on the horizontal trajectories through the phase dia- 
gram, so that 6 = 0.095 in the simulations is mapped to 
S = 0.095 X 0.087/0.137 = 0.0612. 
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